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Abstract 

A robust, user-friendly, and automated method to determine quantum conduc- 
tance in quasi-one-dimensional systems is presented. The scheme relies upon an 
initial density-functional theory calculation in a specific geometry after which 
the ground-state eigenfunctions are transformed to a maximally-localised Wan- 
nier function (MLWF) basis. In this basis, our novel algorithms manipulate and 
partition the Hamiltonian for the calculation of coherent electronic transport 
properties within the Landauer-Buttikcr formalism. Furthermore, we describe 
how short-ranged Hamiltonians in the MLWF basis can be combined to build 
model Hamiltonians of large (>10,000 atom) disordered systems without loss 
of accuracy. These automated algorithms have been implemented in the Wan- 
nier90 code [1] , which is interfaced to a number of electronic structure codes such 
as Quantum-ESPRESSO, Ablnit, Wien2k, SIESTA and FLEUR. We apply our 
methods to an Al atomic chain with a Na defect, an axially heterostructured 
Si/Ge nanowire and to a spin-polarised defect on a zigzag graphene nanoribbon. 
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Wannier function, Wannier90 
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1. Introduction 

Nanostructured materials, such as carbon nanotubes and silicon nanowires, 
promise advances in wide-ranging device applications such as photonics [2], ther- 
moelectrics [3, 4] and biological/chemical sensing[5]. Successful incorporation 
of such structures in real devices requires bottom-up approaches to design, 
which in turn, require an understanding of electronic transport at the nano 
and mesoscales. 
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First-principles simulations based on density-functional theory (DFT) are 
now well-established as a powerful tool for materials modelling. Their success 
is largely due to the high accuracy and computational efficiency that can be 
obtained for many classes of materials. 

Notwithstanding concerns regarding its ability to describe charge transport 
in certain situations [6], DFT combined with the Landauer formulation [7] has 
become a standard starting point for evaluating quantum conductance (QC) [8— 
17]. Calculations typically adopt a 'lead-conductor-lead' geometry (Fig. 1 (top)) 
whereby the conductor is sandwiched between two contacts, or leads, whose 
semi-infinite nature is accounted for by means of surface Green's functions and 
self-energies [18] obtained from standard DFT calculations. 

Despite the success of this approach, realistic nanoscale systems, which typ- 
ically contain arbitrary distributions of impurities, functionalizations and mod- 
ulations of structure and composition are challenging to describe accurately due 
to the asymptotic cubic scaling of conventional DFT calculations with respect 
to system size. 

In this Article, following Lee et al. [14] and Cantele et al. [19], we use a 
method based on the transferability of maximally-localised Wannier functions 
(MLWFs) [20, 21] in order to overcome the cubic-scaling bottleneck. The nov- 
elty of our work lies in the development of robust algorithms for the complete 
automation of the often painstaking manipulations required for preparing a 
Hamiltonian matrix in the MLWF basis. As a result, high-throughput compu- 
tations of QC requiring little user intervention become feasible for disordered 
nanoscale systems. Two further important features of our method are (i) that 
the MLWF basis is optimally compact, ensuring highly efficient determination 
of QC and density of states (DoS) , and (ii) that the nearsightedness of the elec- 
tronic interactions can be exploited in the MLWF basis by piecing together, 
without loss of accuracy, Hamiltonians from DFT calculations on small frag- 
ments to form model Hamiltonians of complex nanostructures consisting of tens 
of thousands of atoms or more. 

The remainder of this paper is structured as follows: Sec. 2 describes briefly 
the underlying theory of Landauer transport and MLWFs, the real-space basis 
in which the transport calculations are performed; Sec. 3 describes the details 
of the implementation of our automated method within the Wannier90 code [1] ; 
in Sec. 4 we present the results of our approach on a number of systems; finally, 
Sec. 5 is reserved for our concluding remarks. 

2. Theoretical Background 

2.1. Landauer Transport 

Within the Landauer formalism, it is assumed that there are no dissipative 
scattering events on the length scale of the conductor region, such that trans- 
mission is coherent, or ballistic. For a single conducting channel at each energy 
E, Landauer showed [7] that the zero-bias, zero-temperature conductance G(E) 
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Figure 1: Top: Schematic illustration of the lead-conductor- lead geometry. Bottom: An 
illustration of how the leads are split into principal layers with Hamiltonian sub-matrices 
labelled according to Eq. (3). 



is given by 



G{E) = ^-T(E), 



(1) 



where T(E) is the probability of transmission through the conducting channel. 
In this framework, G(E) is called the quantum conductance (QC). Extending 
this formalism to multiple channels [22-24] one may write 



G(E) = ^-Tr(T L G r c T R G a c ), 



(2) 



where G^'^ are the retarded (r) and advanced (a) Green's functions associated 
with the conductor, and ^{l.r} are functions that describe the coupling of the 
conductor to the left (L) and right (R) leads. 

The standard approach used to determine the QC of nanostructures that 
has emerged in recent years employs a localised basis set so that the Hamil- 
tonian H of a system in the lead-conductor-lead geometry (Fig. 1 (top)) may 
be partitioned unambiguously. A principal layer [25, 26] (PL) is introduced, 
which is long enough so that (£"|fr|£j n ) ~ if |m — n\ > 2, where Cf is the i th 
basis function in the n th PL and H is the Hamiltonian operator for the entire 
system. By imposing the equality on the Hamiltonian elements, a truncation 
error is introduced, which is controlled systematically by increasing the size of 
the PL (as will be shown in Sec. 4). The Hamiltonian matrix, with reference to 
the bottom panel of Fig. 1, then takes tri-block diagonal form, 
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where interactions between the first PL of the left or right lead and the conductor 
are Hlc and Hcr, respectively; Hf and Hf are matrices formed by orbitals in 
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the same PL in the semi-infinite left and right leads, respectively, and H^ and 
H^ 1 are matrices formed by orbitals in adjacent PLs in the left and right leads, 
respectively. As shown in Eq. (3), these latter four matrices are periodically 
repeated to form Hl and H R (defined in the bottom panel of Fig. 1). 

2.2. Green's Function Formalism 

Knowledge of the seven finite Hamiltonian sub-matrices H £ , H^ a , h^c He 
hen, H R ° and is sufficient to describe the open system of Fig. 1 and to 
calculate the QC from Eq. (2). Following Nardelli [27], in order to determine 
Gq'^ and T{ L R y, we first consider the Green's function G of the whole system, 

(e - H)G = I, (4) 

where e = E + in (77 — ^ 0) for G r . Since G a = (G r )t , in the following we focus 
on G r only and suppress the superscript. From Eq. (4) it can be shown that [18] 

G c = (e-ffc-S L -S fi )- 1 , (5) 

where Sl = h\ c G ®h,LC ancl ^R. = ^ CR G° R hcR represent self-energy terms 
due to the coupling of the conductor to the leads. G°" R y are known as surface 
Green's functions and can be computed efficiently via the iterative procedure of 
Lopcz-Sancho ct al. [28]. Gc is related to the local density of states (DoS) Mc 
of the conductor by [18] 

Mc{E) = --lm(Tv[G c (E)}). (6) 

7T 

Finally, the coupling functions Tr LR \ are given by [18] 

r {i,ii} = *I S {l,_r} - S {l„r}]- ( 7 ) 

In the special case that the lead and conductor are identical and the entire 
lead-conductor-lead system is translationally invariant, the following simplifica- 
tions can be made: Hf = H° R " = H c , and Hf = h LC = h CR = H$. Such 
systems are hereafter referred to as bulk, or pristine, systems and transport 
calculations thereupon are referred to as bulk, or pristine, transport calcula- 
tions. In our results, we will compare the QC of disordered conductors with the 
corresponding bulk, or pristine, QC. 

2.3. MLWF Basis 

In a periodic crystal, within the independent particle approximation, elec- 
trons are described by bands, or Bloch states [V'nk) with band index n and crys- 
tal momentum k. An entirely equivalent representation may be constructed in 
terms of Wannier functions |w n R,), obtained by Fourier transforming |"0 n k) in 
the pair of conjugate variables k and R, where R labels the lattice vector of the 
real-space cell in which the Wannier function is centered. Unlike Bloch func- 
tions, however, Wannier functions may be constructed that are localized in real 
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space, exhibiting exponential decay in systems with an electronic band-gap [29] . 
Even in a metal, exponential localisation can be achieved if an appropriate com- 
bination of filled and empty states is used [14]. 

For an isolated band, observables are invariant under a gauge transformation 
of the form Vnk — > e"^ nk i/> n k. Different choices of the phase <j> n k, however, will 
result in different Wannier functions, and can therefore be used as a means of 
making the resulting Wannier function as localised as possible. Marzari and 
Vanderbilt [20] showed that for a composite, yet isolated group of bands (such 
as those found in the valence manifold of an insulator of semiconductor), one 
may define a set of generalised Wannier functions 



where Umn is a unitary matrix that may be chosen such that the Wannier 
functions are maximally-localised i.e. that the sum of their quadratic spreads 



where (r 2 }„ = (w n0 \r 2 \w n0 } and (r}„ = (w n0 \r n \w n0 ), takes the smallest value 
possible. 

When a set of bands is not isolated from the rest of the band structure by a 
gap across the Brillouin zone, the bands are said to be connected or entangled. 
This is the case in metals and in conduction manifolds of semiconductors and 
insulators. In such cases, within a given energy window, the number of bands at 
each point in k-space varies and the disentanglement procedure of Souza et al. 
[21] is used in order to extract, or disentangle, an optimally-connected subspace 
of a given, constant dimension at each k. Once this optimal subspace has been 
obtained, the usual localization procedure of Marzari and Vanderbilt [20] may 
be applied in order to determine the MLWFs. Once obtained, these provide 
a real-space and often intuitive picture of bonding in materials, to the point 
that they are now used widely as a post-processing tool in electronic structure 
calculations [1]. 

There are a number of advantages to using MLWFs. First, they span a 
much smaller subspace compared to, say, the plane- wave basis in which the 
original ground-state electronic structure calculation is performed. The space 
of MLWFs is arguably the most compact, minimal manifold possible (1 MLWF 
per every band that needs to be described), while still preserving in full the 
accuracy of the electronic structure calculation. As a result, matrices in an 
MLWF basis can be orders of magnitude smaller in each dimension than in 
the original basis, while still reproducing exactly the properties of the ground- 
state, thus enabling very efficient and accurate computation of ground-state 
properties, such as interpolated band structures [30]. For example, the band 
structure of the valence manifold for silicon is equivalently described by ^3000 
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Figure 2: Schematic illustration of the SSG (the single DFT supercell required for automated 
QC calculations). The conductor under investigation is flanked on each side by principal layers 
PL1, PL2 of the leads and a buffer region Bl, B2. The buffer is a length of lead at least as 
large as a principal layer whose function is to ensure the disorder of the conductor has no 
significant effect on the periodicity of the lead Hamiltonian in PL1 and PL2. Also shown is 
the periodic image of PL2 and the regions where each Hamiltonian sub- matrix is derived from 
when expressed in the MLWF basis. 



plane-waves per atom or 2 MLWFs per atom. Second, since they are localised 
in real-space, MLWFs may be used to represent the Hamiltonian of a system 
in sparse matrix form. Finally, this sparsity may be exploited in order to build 
model Hamiltonians of large, structurally complex systems from Hamiltonians 
of smaller fragments. 

3. The Single Supercell Geometry 

The translational symmetry present in crystals is exploited in electronic 
structure calculations by using supercells and periodic boundary conditions 
(PBC). A natural basis set to use for such calculations is that of plane- waves 
and their benefits for DFT calculations are well- understood [31]. 

The lead-conductor-lead geometry of Fig. 1, however, is inherently non- 
periodic. Therefore, as outlined in Sec. 2, if we are to use PBC for our Landauer 
conductance calculations, a transformation to a localised basis set, such as ML- 
WFs, becomes essential. 

Furthermore, the vast range of structural combinations that one could in- 
vestigate means that the change of basis must be coupled to a robust and user- 
friendly algorithm that automatically prepares the Hamiltonian obtained from 
a calculation on a periodic system for use in the transport calculation so that 
high-throughput calculations are possible. The novelty of our work lies in the 
automation of these non-trivial manipulations of Hamiltonian matrices and in 
streamlining the calculations such that a calculation on only a single supercell 
is required; we call this the Single Supercell Geometry (SSG). 

The SSG is shown in Fig. 2, whereby a central conductor is sandwiched be- 
tween a length of lead on the left and right. The conductor is the disordered 
region under investigation and the leads are the contacts whose bulk is peri- 
odically repeated ad infinitum in the open (lead-conductor-lead) system. We 
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split each lead into two parts: the outermost regions must be a PL of lead (PL1 
and PL2) and the inner regions a buffer (Bl and B2) such that any disorder 
within the electronic structure associated with the conductor is localised within 
the region marked Hq. In this respect it is important to converge results with 
respect to the size of the buffer regions; we also impose that BI and B2 must 
be at least one PL of lead in length. 

By transforming to a MLWF basis, our algorithm uses the SSG to identify 
the Hamiltonian sub-matrices required for the transport calculation. Fig. 2 
depicts these regions. For the sake of clarity, it is worth highlighting first that 
Hq is in fact built from the contribution of MLWFs within the conductor and 
the buffers, and second that the interaction between two adjacent PLs of lead, 
H]^ and H^, are built from Hamiltonian matrix elements between MLWFs in 
PL1 and the periodic image of PL2. For this reason, we demand that the left 
and right leads of the SSG be identical in nature. 

The Hamiltonian sub-matrices attained from partitioning the total Hamilto- 
nian require a number of operations performed on them before they can be input 
into transport calculations. First, we need to re-order the MLWFs in real-space 
so that every unit cell in PL1, PL2, Bl and B2 has a consistent sequence of ML- 
WFs. This is because the Hamiltonian corresponding to the semi-infinite leads 
is constructed from sub-matrices extracted from the SSG Hamiltonian in the 
MLWF basis. The connection matrix is constructed from the Hamiltonian 
matrix elements between MLWFs in PL1 and the periodic image of PL2, whereas 
£f£° is constructed from PL1 only. These two matrices are then duplicated along 
the block off-diagonal and block diagonal, respectively, of the Hamiltonian of 
Eq. (3). In doing so, the implicit assumption is that the sequence of MLWFs 
in the rows and columns of the Hamiltonian sub-blocks are the same, which in 
general is not true. To overcome this problem we use the positions of the MLWF 
centres in real-space to order the elements of the Hamiltonian sub-matrices: the 
MLWFs in each unit cell of lead are arranged first according to their position 
along one direction perpendicular to the transport direction, then in the other 
perpendicular direction, and finally along the transport direction itself. This 
ensures that the sub-matrices can be used consistently to build the Hamiltonian 
of Eq. (3). 

The shape of MLWFs are often chemically intuitive and display atomic-like 
or bonding/anti-bonding orbitals. Thus if more than one MLWF exists with 
precisely the same centre, as may happen with rf-like MLWFs on a transition 
metal site, a second level of ordering based on the orbital character of the 
MLWF is performed, employing a technique we have developed using spatially- 
dependent integrals to deduce a unique signature for each MLWF (see Appendix 
A). 

In addition to the ordering of the MLWFs, a second consistency criterion 
must also be imposed. The issue stems from the fact that, although MLWFs 
are always found to be real, they remain undetermined upto an overall sign, 
or parity. As with the issue of ordering the MLWFs, the procedure of building 
the Hamiltonian from sub-matrices implicitly assumes that the MLWFs in PL2 
have the same parity pattern as those in PL1, which in general is not true. 
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To address this issue, we enforce a consistent parity pattern at the level 
of the unit cell of lead onto the ordered Hamiltonian sub-matrices (PL1, Bl, 
B2 and PL2). The parities of the MLWFs in the leftmost unit cell of lead in 
the SSG supercell are used as the template. By assessing the relative parity of 
MLWFs in this unit cell compared to translationally equivalent MLWFs in the 
other unit cells of the PLs and buffer regions, the pattern is enforced throughout 
by multiplying by ±1, as appropriate. The relative parities are determined by 
using the unique signature associated with each MLWF. 

We outline three caveats that apply to the current implementation of the 
SSG method: (i) the Bloch states used as input for determining the MLWF 
basis in the SSG are calculated at the T-point only; (ii) the lattice vectors of 
the SSG must form a orthorhombic set and the direction of conduction must be 
in the x, y or z direction, (iii) the system under investigation must be quasi-one- 
dimensional (although the extension to bulk leads would be relatively simple to 
implement.) 

3.1. Calculation Procedure 

We now outline our general method, this is shown schematically in Fig. 3. 
First we must determine the number of unit cells that make up a PL. Con- 
sider a supercell of lead with 2n + 1 unit cells along the conduction direction, 
whose Hamiltonian in the MLWF basis is found from a T-point DFT calcula- 
tion in PBC. The value of n is chosen such that Hamiltonian matrix elements 
between MLWFs in the central unit cell and the left-most unit cell are less than 
a certain threshold, which is usually set to be around 10 meV. In practice, for 
computational efficiency, instead of a L-point supercell calculation, we apply 
Bloch's theorem to reduce the supercell to a single unit cell, and perform the 
DFT calculation on a regular grid of k-points in the conduction direction. This 
calculation is also used to calculate the bulk, or pristine, QC which is used for 
validation purposes (see Sec. 4). 

Next, the extent of the buffer is determined by assessing the convergence 
of electronic structure in PL1 and PL2 with respect to its size. If the disorder 
present in the conductor region is short-ranged, then often a single PL of lead 
in Bl and B2 is sufficient. This point is illustrated in Fig. 4: we see that beyond 
a few unit cells from a defect (hydrogen functionalization in a (3,3) carbon 
nanotube), the on-site Hamiltonian matrix elements of the MLWFs recover their 
bulk value. Once the extent of the buffer and PLs have been determined, the 
SSG supercell may be built and its Bloch eigenstates found from a conventional 
DFT calculation. This calculation is usually performed in two steps. First, a 
self-consistent calculation at enough k-points to converge the charge density, 
followed by a non-self-consistent calculation at the T-point only, using the self- 
consistent charge density as an input. 

Transformation to the MLWF basis, Hamiltonian-matrix preparation, and 
transport calculations may then be performed. The automated algorithms de- 
scribed above, which are implemented in the Wannier90 [1] code arc designed so 
these steps are performed sequentially, with little or no intermediary user input. 
A natural validation of the results may be performed whereby the disordered 



8 



Automated 
algorithms 
detailed in 
this work 



Preliminary Calculations: 

PL and Buffer Size 



Construct SSG 



PW-DFT Calculation: 



Ground-state density 
Bloch states 



Transform to MLWF Basis: 

(Disentangle if required) 

Real-space Hamiltonian matrix 



Order Hamiltonian 
in MLWF basis 



Enforce Parity Pattern 
onto MLWF Hamiltonian 



Partition MLWF Hamiltonian: 
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Figure 3: A flow diagram depicting the key steps in our calculation procedure. 
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Figure 4: Illustration of the inherent electronic nearsightedness in a (3,3) carbon nanotube 
functionalised with a single hydrogen atom. The white-gray spheres represent the atoms of 
the structure while the colored spheres represent the deviations from the "bulk" values of the 
on-site Hamiltonian matrix elements for each MLWF. The size of colored spheres is another 
indication of the deviation of the matrix element from its "bulk" value. The smaller the 
sphere, the smaller the deviation. 

conductor region of the SSG is replaced by a section of pristine lead: identical 
results should be achieved with the bulk calculation. 

3.2. Combination of multiple defects 

Moving from a Bloch to a Wannier representation is not only a means by 
which to represent electronic structure in a very compact manner. It also opens 
the possibility to exploit the real-space nature of the basis to build very large 
systems - systems so large that a conventional DFT calculation would be in- 
tractable. 

The fact that electronic nearsightedness becomes explicitly manifest in the 
MLWF basis, as highlighted in Fig. 4, allows them to be used to build the 
Hamiltonian matrix of a large structure from the smaller Hamiltonian matrices 
of its constitutive sub-systems. 

In order to illustrate the method, consider the schematic lead-conductor-lead 
system shown in Fig. 5 in which the conductor region has two identical defects 
separated by a region of lead material in the form of a buffer (Bl' and B2'). 
We could calculate the QC of this structure by making a SSG with the whole 
conductor (regions X and Y). However, we may exploit the nearsightedness 
of the MLWF basis to find a more computationally efficient approach. If the 
effect of the defects is localised (in the sense that the local electronic structure 
and geometry at the junction between Bl' and B2' is sufficiently similar to 
that seen in the leads), then we may construct the Hamiltonian for the system 
with two defects (Fig. 5) from information gathered from one SSG calculation 
containing just a single defect (Fig. 6). Since this system is smaller, there is a 
clear advantage in terms of computational cost for the initial DFT calculation. 

The Hamiltonian of the conductor in the two-defect system may be written 

as 



where X, Y and XY represent blocks of Hamiltonian matrix elements among 
MLWFs in region X, among MLWFs in region Y, and between MLWFs in these 
two regions, respectively. 
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Figure 5: Schematic for a SSG with a conductor containing two identical defects. We identify 
two regions in the conductor, X and Y. 
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Figure 6: Schematic for a SSG with a conductor Z containing one defect. 



Given the geometry of the system, and the nearsightness of the electronic 
structure, blocks Hx and Hy should be quite close in terms of their matrix 
elements. Moreover, because of the constraint that a buffer is at least as large 
as a lead principal layer, we expect the non-zero matrix elements of Hxy to 
correspond closely to the overlaps between the two adjacent principal layers. 
This observation stems from the very definition of a principal layer. As a con- 
sequence, we can construct a close approximation to He by using the matrices 
extracted from a SSG calculation of the structure shown in Fig. 6. In this ap- 
proximation, blocks Hx and Hy are replaced with Hz-, and Hxy is replaced by 
the overlap matrix between two principal layers of lead (namely Hj^ , see Fig. 
2), i.e., 

%\ (") 

An example of this approach is demonstrated in Sec. 4 for a defected silicon 
nanowire. 

The approach described above is general and may be applied to any number 
of isolated defects in the conductor region. In this way, Hamiltonians for systems 
of almost arbitrary size may be constructed with first-principles accuracy from 
one DFT calculation in a SSG with a single defect. We note in passing the 
importance that the MLWFs parities are consistent between different regions of 
the system. As mentioned in Sec. 3, the parities need to be checked and made 
consistent to allow seamless connections between Hamiltonian sub-matrices, a 
task that is automatic in the present approach. 

Furthermore, the Hamiltonian of a conductor with more than one type of 
defect may be constructed by combining matrix elements from separate SSG 
calculations. In this latter case, care must be taken in order to align the Fermi 
energies of the two (or more) distinct calculations. This is the consequence of 
the lack of an absolute reference for the electrostatics in PBCs, which can lead 
to Fermi energies that are shifted by a constant. 
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Additionally, building a large-scale structure with tens of defects and tens of 
thousands of atoms can be a painstaking task. In order to simplify this process, 
we have designed a utility package to the Wannier90 code that helps the user 
create these large scale structures. From a single Wannier90 calculation in the 
SSG geometry, both randomised and custom-made structures can be built with 
case. An illustration of the use of that functionality is given in the fourth 
example of Sec. 4. 

4. Applications 

We present now a number of examples using the method described in Sec. 2. 
The aim is to illustrate its robustness in a range of applications: beginning 
with a defected atomic chain, and building complexity via a heterostructured 
nanowire and a spin-polarised graphene nanoribbon. Finally, we provide an 
example to validate the use of SSG Hamiltonian fragments in the construction 
of model Hamiltonians for larger systems. All DFT calculations are performed 
with the Quantum-ESPRESSO package [32] and with (unless otherwise stated) 
norm- conserving pseudopotentials. [33] 

Atomic Al Chain 

First, we consider the QC of an Al chain with a single Na atom substitutional 
defect. The construction of the SSG is performed with care: a suitable PL 
length must first be decided upon by assessing the rate of the decay of the 
matrix elements of the Hamiltonian between MLWFs. Additionally, the defect 
is expected to have a large effect on the electronic structure, thus the buffer size 
must also be carefully chosen. 

To assess the length of a PL we use the method outlined in the Sec. 3.1, 
whereby the Hamiltonian in the MLWF basis of a single unit cell of lead is de- 
termined at many k-points. We perform the DFT calculation on a single unit 
cell (consisting of one Al atom) , with a regular grid of 32 k-points along the ex- 
tended direction. The total energy is converged to 10 -11 eV using a 500 eV en- 
ergy cut-off, exchange and correlation are described by the PBE functional [34]. 
The unit cell is 2.47 A long in the conduction direction, with 10 A separating 
periodic images in the transverse directions. We proceed to the determination 
of the MLWF basis by disentangling three Wannier functions from 30 bands. 

Fig. 7 shows the decay of the interaction between MLWFs by averaging 
the on-site Hamiltonian elements (w n o\H\w n ji) between equivalent MLWFs in 
different unit cells labelled by the primitive lattice vector R (black solid line). 
The maximum matrix element (red crosses) gives the maximum error incurred 
due to truncation of the interaction if the PL were to be cut at that unit cell. In 
this case, the PL is chosen to be eight unit cells long, such that the maximum 
truncation error is approximately 11 meV and the average error is 9 meV. 

A buffer size of one PL plus three unit cells is chosen for the SSG such that 
on relaxation the RMS difference in the position of the MLWF centres from their 
ideal bulk position in the rightmost unit cell of PL1 is less than 5 x 10 -3 A. 
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Figure 7: Decay of the Hamiltonian elements {w n o\H\w n n) between increasingly distant Al 
unit cells in units of the primitive lattice vector. Cell averaged elements are shown in black 
(error bars show a standard deviation); the largest Hamiltonian values between unit cells 
increasingly distant from R. = are shown by red crosses. The dashed line highlights the 
chosen PL size (see text). 

The SSG therefore consists of a total of 39 atoms. Performing the transport 
calculation provides the QC and DoS shown in Fig. 8 (red, dashed; centre and 
right panels, respectively). For comparison, the bulk band structure, QC and 
DoS (black; left, centre and right panels, respectively) are also shown. In the 
bulk case there are clear contributions from the s band and two degenerate p 
bands to the QC: these are both significantly reduced in the defected case, with 
conductance close to zero at lower energies. This may be interpreted as the 
hybridization of the s orbital associated with the Na to adjacent Al p orbitals. 

4-2. Si/Ge Nanowire Helero structures 

We now increase the complexity of the SSG system by considering a thin 
(0.39 nm radius) Si nanowire in the (110) direction with a Ge heterostructure 
inserted as a defect. The DFT calculations detailed in this example (and those 
on the nanowires of in Sec. 4.4) are performed within the LDA and with an 
energy cutoff of 400 eV. We begin with a single cell for the lead (8 Si atoms, 8 
H atoms) (see Fig. 9 (top)) and perform a DFT calculation with 20 k-points, 
allowing atomic positions and lattice parameter in the conduction direction to 
relax. Forces are converged to 5 meV/A. Once the ground-state is found, we 
transform to the MLWF basis and assess the PL length in the usual manner. 
Fig. 10 displays the decay of the Hamiltonian matrix elements as a function of 
increasingly distance, using the same notation as Fig. 7. With four unit cells in 
a PL the average truncation error is below 2 meV. 

A SSG is built by repeating the single unit cell of Si and inserting 5 copies 
of a similarly relaxed Ge unit cell (see Fig. 9, bottom panel). Without further 
relaxations of the geometry, it was found that a single PL was sufficient to 
converge the electronic structure in PL1 and PL2 to that of a bulk lead. Hence 
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Figure 8: Left: Band structure of the single cell, single Al atom bulk transport calculation. 
Centre: QC comparison of bulk (black, solid) and defected SSG (red, dashed) systems. Right: 
DoS comparison of bulk and defected SSG. 




Figure 9: Top: Front and side view of single Si cell used for PL determination and bulk 
transport calculations. Bottom: Supercell of the Si/Ge nanowire for use in the SSG method. 
Red, cyan and magenta represent H, Si and Ge atoms respectively. 

our SSG consisted of 16 Si unit cells with five Ge unit cells sandwiched at their 
centre (Fig. 9, bottom panel). Using our automated routines, this 336 atom 
unit cell provides the valence QC and DoS shown in Fig. 11 (red, dashed lines; 
centre and right panels, respectively). The bandstructure, QC and DoS of the 
pristine silicon nanowire are also shown (black solid lines; left, centre and right 
panels, respectively). The drop-off in conductance just below the Fermi level 
is due to localization of the highest occupied molecular orbital within the Ge 
quantum well. 

4-3. Spin-polarised graphene nanoribbon 

In this third example, we look at a spin-polarised graphene nanoribbon func- 
tionalised with a single hydrogen atom. As with the previous examples, we start 
with a calculation on a single unit cell. Our system is a zigzag nanoribbon of 
length 2.46 A with a width of 9.27 A. A regular grid of 20 k-points is used in 
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Figure 10: Decay of the Hamiltonian elements {w n o \H\w n n) between increasingly distant 
Si nanowire unit cells. The notation used is equivalent to that of Fig. 7, where average 
on-site elements are shown in black (error bars show a standard deviation) and the largest 
Hamiltonian value between unit cells increasingly distant from unit cell are shown by red 
crosses. The dashed line highlights the chosen PL size (see text). 
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Figure 11: QC (centre panel) and DoS (right panel) for pristine Si nanowire (solid, black 
lines) and axially heterostructured Si/Ge nanowire (dashed, red lines). The bandstructure of 
the pristine silicon nanowire is also shown (left panel). 
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Figure 12: Interpolated band structure (up spins) for a pristine zigzag graphcnc nanoribbon 
(black solid lines), and "exact" band structure in the complete plane- wave basis set given by 
the electronic structure code (red dots). The Fermi energy is set to eV. 

the conduction direction and the supercell is built such that a vacuum region 
of 10 A lies between periodic images. A cutoff of 400 eV for the kinetic energy 
and 4500 eV for the charge density is used together with a PBE exchange and 
correlation functional and ultrasoft pseudopotentials[35]. Both the atomic posi- 
tions and the unit cell length were fully relaxed; individual forces are less than 
18 meV/A. 

We must specify a starting non-zero magnetization such that the self-consistent 
loop converges to a magnetic state (in this case we restrict ourselves to a ferro- 
magnetic state, even though the ground state is anti- ferromagnetic [36]). Next, 
a non-self consistent calculation is used to compute the band energies for both 
"up" and "down" spins. It is important at this stage to compute a sufficient 
number of bands to capture the entire ir manifold, otherwise the disentangle- 
ment of the conduction manifold would be meaningless. In our calculations, we 
used 30 bands and we kept all the bands up to —0.5 eV in the frozen window 
for the disentanglement procedure. 

The quality of the disentanglement may be assessed by comparing the inter- 
polated band structure provided by Wannicr90 to the full band structure given 
by the electronic structure code (for this ferromagnetic state, both spin types 
have an almost identical band structure except around the Fermi level). As can 
be seen from Fig. 12, the match between the Wannier interpolation (solid lines) 
and the "true" band structure (dots) is excellent. We see that the interpolated 
band structure describes perfectly conduction states upto about 2.5 eV above 
the Fermi energy. 

Satisfied with the MLWF transformation, the PL size is assessed in the 
manner described earlier. Choosing a PL size of four unit cells (with a maximum 
truncation error of less than 46 meV) is sufficient for this example. 

The SSG consists of 4 unit cells of lead in both the PL1, PL2, Bl and B2 
regions. Two unit cells form the conductor region (see Fig. 13), upon which a 
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Figure 13: Top and side view of the SSG structure used to compute the quantum conductance 
of both "spin up" and "spin down" channels. Carbon atoms are yellow and hydrogen is shown 
in red. The overall supercell consists of 18 unit cells of lead and an extra hydrogen, which 
totals 181 atoms. 

hydrogen atom is placed at an approximate C-H bond length above one of the 
carbon atoms. The whole system is fully relaxed both for atomic positions and 
supercell length in the direction of conduction. The force convergence criteria, 
cutoffs and energy convergence criteria are similar to the ones used in the single 
lead unit cell case above. 

After relaxing the structure, the next step is to perform a spin-polarised DFT 
calculation followed by a non-self consistent calculation to extract the "up" and 
"down" band energies at T. The QC is calculated in the usual manner for each 
spin channel separately with Wannier90. The result for the spin-dependent QC 
is shown in Fig. 14. 

We see on the graph that at the Fermi level the quantum conductance of 
the system is slightly spin-polarised with a majority of "up" spins. Applying a 
slightly negative bias, we see that the system can be close to 100% spin-polarised. 
The opposite spin polarization can be achieved with a slightly positive bias. 

4-4- Doubly Defected Si Nanowire 

This final example demonstrates an extension to the SSG method whereby 
the sub-Hamiltonians it creates are manipulated and combined to construct 
model Hamiltonians of larger systems (see Sec. 3.2). All DFT calculations out- 
lined here are performed within the LDA, at the T point, with the same energy 
cut off and tolerances described in Sec. 4.2. 

The system we aim to describe is shown in Fig. 15 (top): a Si nanowire 
with two single cells of Ge separated by a length of Si. This system may be 
thought of in two ways: first, as a Si nanowire with a single defect containing 
the two cells of Ge and the separating Si cells; and second as doubly defected 
Si nanowire, with each defect being a single cell of Ge. These two perspectives 
lead to two methods by which we can determine the QC. The first suggests a 
SSG calculation in which the conductor region contains the two Ge defects. The 
resulting QC is see in solid black in Fig. 15 (bottom). This is compared to the 
QC derived from a calculation in which the multiple defect method described 
in Sec. 3.2 is used. 

For the doubly defected case, we perform a SSG calculation on a Si nanowire 
with a single Ge cell defect and use the Hamiltonians provided to build a Hamil- 
tonian of the larger system in question (Fig. 15 (top)). The QC from this 
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Figure 14: Spin-dependant QC close to the Fermi level. One can clearly see that depending 
on the applied bias, a spin-polarised current can be induced in this system. 

calculation is seen in dashed red in Fig. 15 (bottom). The two calculations 
agree remarkably well, validating the premise of nearsightedness discussed ear- 
lier. Since the computational expense of conventional DFT methods scales as 
0(N 3 ), where N is the number of atoms in the supercell, the multiple de- 
fect method represents a significant step forward to describe realistic system 
sizes with first-principles accuracy, and can be used to construct faithful model 
Hamiltonians for systems that contain tens of thousands of atoms[37, 38]. 

5. Conclusions 

In this paper we have presented a user-friendly and automated approach to 
calculate the quantum conductance and density of states in quasi-one-dimensional 
systems. The method converts the Bloch eigenstates of a single DFT calcula- 
tion, within our single supercell geometry, to a basis of MLWFs. In this basis we 
determine the electronic transport properties by automatically extracting the 
Hamiltonian sub-matrices required for the transport calculation. To illustrate 
the robustness, wide applicability and efficiency of our method, we have pre- 
sented calculations on an atomic Al wire, a spin-polarised graphene nanoribbon, 
and axially heterostructured Si/Ge nanowires. Furthermore, we have shown how 
the transport properties of meso-scale conductors that are beyond the current 
capabilities of conventional first-principles electronic structure calculations can 
be calculated with first-principles accuracy by exploiting the transferability of 
MLWFs as building blocks of large model Hamiltonians. 
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Figure 15: Top: Si nanowire with two Ge heterostructure defects. The system is investigated 
by manipulating the Hamiltonian of a single defect (multiple defect method) and directly using 
a SSG. Red, cyan and magenta atoms are H, Si and Ge respectively. Bottom: Comparison of 
QC for the two methods, showing excellent agreement. 
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Appendix A. MLWF Signatures 

Here we detail the set of spatially-dependent integrals that we use to deter- 
mine a signature for each MLWF. These signatures are used for two purposes. 
First, they enable a sorting algorithm to distinguish between MLWFs of differ- 
ent shapes with similar centers. Thus they may be ordered consistently over 
between unit cells - a key requirement for our approach. Secondly, they are used 
to determine the relative parity of MLWFs so that a consistent parity-pattern 
may also be enforced. 

We begin with the integral 




(A.l) 
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where V is the volume of the cell, q is a vector in reciprocal space and r c is 
the centre of Wannier function w n (r) (we assume sampling at T-point only). 
One may write jfl n (r) = J2 m U m nU m {r), where u m (r) is the periodic part of the 
Bloch wavefunction at band m. U mn is the unitary matrix found in Eq. (8) that 
minimises the spread of the Wannier functions. u m (r) can be written in terms 
of its discrete Fourier transform u m (G), u m (r) — ^ G u r „(G)e 4G ' r . Thus, the 
integral in Eq. (A.l) may be written as 



Jn(q) 



(A.2) 



where q is a G-vector of the form lh\ + nib? + nh^, where {/, m, n} G Z and 
{bi,b2,b3} are the reciprocal lattice vectors. Equating real and imaginary 
parts of Eq. (A.l) and Eq. (A.2), one may write 



C(q) = 



^ J w n (r) cos(q • (r - r c )) dr 

' :q - rc ]Tt/ m „«; n (q) , 



v 

Re 



(A.3) 



and 



C(q) 



v 



Im 



w n (r) sin(q- (r - r c )) dr 



(A.4) 



Since most DFT codes compute u m (G), obtaining any set of I n incurs negligible 
computational expense. 

The set of integrals that are used to determine a signature are given by 

I n = ^ J w n (r) sin" (jr-(x - x c )^j sin' 3 (j^-(y - Vc)j sin 7 (jr-(z - z c )^j dr 

(A.5) 

where r c = (x c , y c , z c ), V = L x L y L z , a, /3, 7 G {0, 1, 2, 3} and a + (3 + 7 < 3. 
Each of the resulting 20 integrals may be written as linear combinations of 
those outlined in Eqs. (A.3) and (A.4). The signature of the MLWF is thus 
given by the 20-element unit vector of these integrals. Dot products between 
two MLWFs' signatures reveal in a compact form their relative shape and parity. 
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